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Abstract — An easy-to-use and effective formula for stability 
testing of a system witli fractional-delay characteristic equation in 
the general form of A(s) = Po{s) + J2Zi P^i^) exp(-Os'^') = 0, 
where Pi{s) (i = 0, . . . , TV) are the so-called fractional-order 
polynomials and (i and Pi are positive real constants, is proposed 
in this paper. The proposed formula determines the number of 
unstable roots of the characteristic equation (i.e., those located 
in the right half-plane of the first Riemann sheet) by applying 
Rouche's theorem. Numerical simulations are also presented to 
confirm the efficiency of the proposed formula. 

L Introduction 

In some of the recently-developed control problems we need 
to check the stability of a system with the so-called /racf/ona/- 
delay characteristic equation in the general form of 

N 

A{s) = Po{s) + P^{s) exp(-Gs^') = 0, (1) 



where Ci and Pi are positive real constants and Pi{s) {i = 
1, . . . , N) are fractional-order polynomials in the form of 

Mi 

=^a,fcs"'^ (2) 
fc=i 

where Uik and aik are real and positive real constants, respec- 
tively, and 



(3) 



where, without any loss of generality, it is assumed that the 
powers of s in ([3]l satisfy the following relations: 



a„ > > • • • > ai > 0. 



As an example of a system with fractional-delay characteristic 
equation, consider a classical unity-feedback system in which 
a process with transfer function 



Gis) 



1 + sT 



(5) 



is controlled with the so-called PI D** controller with transfer 
function |[T]: 



C{s) =Kp[l + 



1 



(6) 



where Kp, Ti, Td, A, and /i are unknown parameters of the 
controller to be determined. It can be easily verified that the 
characteristic equation of this system is as the following 

A(s) = T,s^{l + sT)+KpK{TiS^ + l + TiTdS^+'')e--'^ = 0, 

(7) 

which can be considered as a special case of ([T]i with iV = 1, 

Pois) = T,s\l + sT), (8) 



(9) 



Ci = L, and /3i = 1. If in this example one tries to 
find the optimal values of Kp, Ti, Td, A, and /i by means 
of meta-heuristic optimization algorithms such that a certain 
cost function (e.g., ISE performance index corresponding to 
the tracking of unit step command) is minimized, he/she 
will need a method to check the feasibility of the solutions 
generated by the meta-heuristic optimization algorithm from 
the stabiUty point of view. It should be noted that since in such 
optimization problems the cost function is usually expressed 
in the frequency domain (by applying Parseval's theorem), 
the resulted optimal controller may destabilize the feedback 
system |2|. 

As a more complicated example, consider the problem of 
designing an optimal PI'^'D^ controller for a process whose 
transfer function consists of fractional powers of s possibly in 
combination with exponentials of fractional powers of s. For 
example, the transfer functions: 



(4) and 



G{s) 



Gis) 



y/s sinh (y/s) ' 
sinh (xoy/s) 



< xo < 1, 



< xo < 1, 



(10) 



(11) 



sinh (\/s) 

appear in boundary control of one-dimensional heat equation 
with Neumann and Dirichlet boundary conditions ||3l. Other 
examples of this type can be found in |3|-|5l. Moreover, in 
some applications in order to arrive at more accurate models, 
the process is modelled with a fractional-order transfer func- 
tion. For instance, Podlubny [SI showed that the fractional- 
order transfer function: 

^('*) ~ n -7^1^0.2.5708 ^ rr,oor_ns^79 , 1 rrr^n^ (1^) 



0.7943s 



5.2385s' 



,0.8372 



1.5560' 



can better model a heating furnace compared to classical 
integer-order transfer functions. Clearly, in dealing with com- 
plicated transfer functions such as those given in (fT0t-(fT2]i 
we need more powerful tools to determine the stability of the 
corresponding closed-loop system. 

Stability analysis of the feedback system when such com- 
plicated transfer functions exist in the loop is a challenging 
task. Even the stability analysis of a feedback system which 
consists of both PI'^D'^ controller and a process with dead- 
time is not straightforward. So far, many researchers have tried 
to develop analytical or numerical methods for stability testing 
of systems with fractional-delay characteristic equations (see 
Q for a detailed review of some important works in relation 
to the stability testing of fractional-delay systems). Probably, 
the most famous analytical method for stability testing of 
fractional-order systems (as a special case of fractional- delay 
systems) is the sector stability test of Matignon fS*!, which 
was akeady reported in the work of Ikeda and Takahashi |j9j| . 
AppUcation of this method is limited to the case where the 
sigma term does not exist in ([T]i and Po(s) is of commensurate 
order. Few numerical algorithms for stability testing of ([T) 
can also be found in the literature (see, for example, IfTOl and 
Q and the references therein for more information on this 
subject). As far as we know, all of these methods suffer from 
the limitation that can be applied only to a certain class of 
fractional-delay systems ifTOl . or the results are of probabilistic 
nature [TJ. 

The aim of this paper is to propose a formula for determin- 
ing the number of unstable roots of ([T). The proposed formula 
is actually a generalization of the method already proposed by 
author in ifTol . However, the formula developed in this paper 
has the advantage of being much simpler compared to the one 
presented in [1P|, and moreover, it can be easily applied to a 
more general form of fractional-delay characteristic equations. 

The rest of this paper is organized as follows. The proposed 
formula for stability testing of fractional-delay characteristic 
equations is presented in Section Four numerical examples 
are studied in Section |III] and finally. Section |IV] concludes 
the paper 

II. Proposed formula for stability testing of 

FRACTIONAL-DELAY CHARACTERISTIC EQUATIONS 

The first step in dealing with multi-valued complex func- 
tions (such as the one presented in ([T]i) is to construct the 
domain of definition of the function appropriately. The domain 
of definition of the characteristic function given in ([TJ is, 
in general, in the form of a Riemann surface with infinite 
number of Riemann sheets, where the origin is a branch point 
and the branch cut is considered (arbitrarily) at Mr . Equation 
A(s) = as defined in ^ has, in general, infinite number of 
roots which are distributed on this Riemann surface. As a well- 
known fact, a system with characteristic equation ([1) is stable 
if and only if it does not have any roots in the right half-plane 
of the first Riemann sheet jTl, ifTTI . Hence, stability analysis 
of a system with characteristic equation ([T]i is equivalent to 



investigation for the roots of A(s) = in the right half- 
plane of the first Riemann sheet. In the following we will 
use Rouche's theorem for this purpose. 

First, let us briefly review the Rouche's theorem. Consider 
the complex function f : C C which has zeros of orders 
nil, . . . , ruk respectively at zi, . . . , Zk and does not have any 
poles. This function can be written as 

f{s) = g{s){s - zi)™i X (s - 22)"^ X • • • X (s - Zfe)"'' , (13) 

where g(s) has neither pole nor zero. Taking the natural 
logarithm from both sides of (fT3] l leads to 

ln/(s) = In 17(5) + TOi ln(s — Zi) + m2 ln(s — Z2) + . . . 

+ mfeln(s-Zfe). (14) 

Derivation with respect to s from both sides of (fT4l i yields 

fjs) ^ g'js) ^ ^ m2 ^ ^ vtik ^^^^ 

f{s) g{s) s- zi s- Z2 "' s- Zk 

Now let 7 be a simple, closed, counterclockwise contour such 
that /(s) has no zeros (or singularities like branch point and 
branch cut in dealing with multi-valued functions) on it. Then 
it is concluded from the Residue theorem that 

where M is equal to the total number of the roots of f{s) = 
inside 7. Clearly, if the contour 7 is considered such that 
all zeros of /(s) lie inside it then we have M = "^j- 
Equation dTSI ) can be used to calculate the number of zeros of 
the given function /(s) inside the desired contour 7 (which, of 
course, should have the above-mentioned properties). For this 
purpose, we can simply use a numerical integration technique 
to evaluate the integral in the right hand side of (fT6t for the 
given contour 7 and function /. 

According to the above discussions, by setting f{s) equal 
to A(s) and 7 equal to the border of the region of instability 
(which is equal to the closed right half-plane of the first 
Riemann sheet) the value obtained for M from ( fTSI l will be 
equal to the number of unstable roots of the characteristic 
equation. In the following, we consider the contour 7 as shown 
in Fig. [T] and f{s) — A(s) (where A(s) is defined in ^) 
and then simplify the integral in the left hand side of (fTSI l 
to arrive at a more effective formula for stability testing of 
the fractional-delay system under consideration (clearly, the 
system is stable if and only if 1\I = 0). Note that the very 
small semicircle in Fig. [T] is used to avoid the branch-point 
located at the origin. 

According to (fTST i and Fig. [T] we can write 




Fig. 1. The contour 7 considered on the first Riemann sheet for stability 
testing of the fractional-delay system under consideration. A system with 
characteristic equation A(s) = (as defined in {!)) is stable if and only if it 
does not have any roots inside 7. 



In ( [TtI i. the integral J^^^^^ is calculated as 

= -/ —^^du;+ ^ . ^ Hd^ (19) 

which yields 

The integral J^^ in (fTTI i is calculated as 

In the above equation lime_i.o A(£e*^) is equal to a nonzero 
constant (else, the characteristic function has a strong singu- 
larity at the origin and the corresponding system is unstable) 
and A'(ee'^) ^ Ke^ as e — > 0, where K and rj > —1 are two 
constants. Hence, f tends to zero as e 0. Finally, f in 
([TtT i is calculated as 

/ = lim r ^^y^ ffie'^dg (24) 

= ia„7r. (26) 

(See ([TJ and ([3j for the definition of a„.) Substitution of (1211 1 
and ( |26l ) in ([IT) and considering the fact that J =0 results 

M.^-ir Re|§i4|.„. (27, 



where Af is equal to the number of unstable poles of a system 
with characteristic equation A(s) = as defined in ([T]). 

Equation dZTI l is the main result of this paper It should be 
noted that the value of e in dZTl i cannot, in general, be consid- 
ered exactly equal to zero. That is because of the fact that the 
numerical integration technique used to evaluate the integral 
in ( I27] ) performs this task by evaluating the integrand at 
different points of the lj axis. Hence, the numerical integration 
algorithm may halt if the integrand becomes singular at the 
origin (which is the case if, for example, < ai < 1 in (O). 
In practice, in order to determine the number of unstable poles 
of the given fractional-delay transfer function we can consider 
the lower and upper bound of the integral in (l27t equal to 
sufficiently small and big positive numbers, respectively. The 
MATLAB function quadi (as well as quadgk) can be used to 
evaluate the integral in (|27T i. Some numerical examples will 
be presented in the next section. 

III. Numerical examples 

In the following we study the application of dZTl ) for 
stability testing of some fractional-delay systems. In each 
case, the impulse response of the corresponding system is 
also plotted to verify the correctness of the result. The 
method used in this paper to calculate the impulse response 
of the given fractional-order system is based on the for- 
mula proposed in |12| for numerical inversion of Laplace 
transforms. In this method the impulse response of the 
given fractional-order system is approximated by numeri- 
cal inversion of its transfer function. The MATLAB code 
of this method, invlap.m, can freely be downloaded from 
http://www.mathworks.com/matlabcentral/fileexchange/ Most 
of the following examples have already been studied by author 
in fll. 

Example 1. Consider a system with characteristic equation 

Ai(s) = + (28) 

= ^STT/e ^ ^7r/2 ^ ^7r/3 _^ ^ = Q. (29) 

The roots of this equation can be calculated analytically, which 
are as the following: 

^ ^M2k^ + l) ^ fciGZ, (30) 

and 

= e'^'-^'''+^\ k2eZ. (31) 

As it is observed, the characteristic equation given in (|29) has 
infinite many roots which are distributed on a Riemann surface 
with infinite number of Riemann sheets. It is concluded from 
( l30b and ( |3TI ) that ( |29] l has four roots on the first Riemann 
sheet which are e^^'^ and e^^^, and none of them are located 
in the right half-plane (recall that all roots whose phase angle 
lies in the range [— tt, tt) belong to the first Riemann sheet). 

Comparing (|29] l with ([TJ and (|3j yields a„ — hir/Q (note 
that ( |29] l has no delay terms). Application of (|27| i assuming 
that the lower and upper bound of integral in (l27t are equal to 
and 1000, respectively, yields M = 2.3300 x IQ-* which is 
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Fig. 2. Impulse response of a system with ti'ansfer function i32\ . 



Fig. 3. Impulse response of a system with transfer function 1 134) . 



consistent with the above-mentioned analytical result. Figure |2] 
shows the impulse response of a system with transfer function 

1 1 



(32) 



Ai(s) 

As it can be observed in this figure, the impulse response of 
the system is absolutely summable, as it is expected. 

Example 2. Stability of a system with fractional-delay 
characteristic equation: 



A2(s) = s + X(^^+l)e- 



(33) 



is studied in [TSl and it is especially shown that it is stable 
for K < 21.51 and unstable for K > 21.51. Application of 
(|27] | assuming that K = 21, — 1, and the lower and upper 
bound of integral are equal to and 500, respectively, yields 
M = 3.4227 X 10^^, which implies the stability of system 
as it is expected. Figure |3] shows the impulse response of a 
system with transfer function 

'''^'^^A^) = s + 2liJ+l)e-^- ^^^^ 

As it can be observed, the impulse response is absolutely 
summable, as it is expected. Repeating the above procedure 
with K ^ 22 yields M = 2.0174, which means that in this 
case the system has two unstable poles. This result is also 
consistent with the one presented in iflOl . 

Example 3. It is shown in [14J that a system with charac- 
teristic equation 



Ms) 



,1-5 



1.5s + 4s' 



0.5 



8- 1.5se- 



0, (35) 



is stable for the values of r e (0.99830, 1.57079) and unstable 
for other values of t. It is also shown by author in |10| that 
this system has two unstable poles for r — 0.99. Application 
of dZTl i assuming t — 1 and considering the fact that here we 
have Q!„ = 1.5 yields M = 0.0082 (the lower and upper bound 



= 0.1 
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Fig. 4. Impulse response of a system with transfer function ([36). 



of integral are considered equal to and 500, respectively). 
As it is observed, the result obtained by using the proposed 
method is fairly close to zero. Figure shows the impulse 
response of a system with transfer function: 

H^is) = — ^ = -pp . (36) 

^' Asis) 1.5s + 4sO-5 + 8- 1.5se-^ 

As it can be observed in this figure, the impulse response 
of the system is absolutely summable and consequently, the 
corresponding system is stable. In this example, application 
of dZTl i assuming r — 0.99 yields 1.9994 which is consistent 
with the result presented in ifTol . 

Example 4. It is shown in |7| (by applying Lambert 
W function) that a system with the following characteristic 
equation 



is stable. Clearly, here we have 



+ = 0, (37) 
5/6. Application of 
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Fig. 5. Impulse response of a system with transfer function (|38). 



dZTl l (assuming that the lower and upper bound of integral are 
equal to and 100, respectively) leads to M = 0.0290, which 
implies the stability of system. Figure |5] shows the impulse 
response of a system with transfer function 



It can be observed that the impulse response is absolutely 
summable and consequently, the system is stable, as it is 
expected. 

IV. Conclusion 

An easy-to-use, effective and very general formula for 
stability testing of fractional-delay systems is proposed in this 
paper. The proposed formula can be used to determine the 
number of unstable poles of a system whose characteristic 
equation contains, in general, both fractional powers of s and 
exponentials of fractional powers of s. 
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